-module(img).

-export([run/0]).

% enable    native compiler
-compile([debug_info, {hipe,o3}]).


-record(gray_vals,{w,h,byte_tuple}).

run() ->
%	ImgFile = "q-man.jpg",
%	ImgOut = "n_edge_avg40_20_10_5_2.png",
%	W=500,
%	H=500,

	ImgFile = "blotch.png",
	ImgOut = "blotch_e2.png",
	W=308,
	H=300,

	%ImgFile = "q-man-gray-div2.jpg",
	%ImgOut = "medgeD2All.png",
	%W=250,
	%H=250,

    %ImgFile = "q-man-gray-div5.jpg",
	%W=100,
	%H=100,
	
    %ImgFile = "roadsign.svg.png",
	%ImgOut = "medgeD5All3Road.png",
	%W=320,
	%H=320,
	
	

	PrivDir = code:priv_dir(imgsegm) ++ "/",
	Grays = list_to_tuple(port_recv_all("convert " ++ PrivDir ++ ImgFile ++ " gray:-")),
	io:format("recv size ~p~n",[size(Grays)]),
	GrayVals = #gray_vals{w=W,h=H,byte_tuple=Grays},
	%log_block(GrayVals, 200, 309, 10),
	%log_block_subpixel(GrayVals, 200, 309),
	%{200,309,radius,10},
	%write_img_lines(GrayVals,309),
	write_img_edges(GrayVals,2),
	%{line_edges(GrayVals,309,2),
	% rotated_line_edges(GrayVals,5,{200,309}, {1,0.05})  },
	%write_img("medge01.png",W,H,  all_rotated_line_edges_gray_fun(GrayVals,2,-0.1)).

	% sharp edges - but many edges are not sharp!
	%ESetAll = rotated_line_edges_coalesce(GrayVals,2,5),
	
	AllRadians = all_radians(5),
	%EDict12 = rotated_line_edges_coalesce(GrayVals,12,AllRadians),
	%EDict8 = rotated_line_edges_coalesce(GrayVals,8,AllRadians),
	%EDict = dict_sum(EDict8, EDict12),
	%EDict = rotated_line_edges_coalesce(GrayVals,40,AllRadians),
	EDict = rotated_line_edges_coalesce_sum(GrayVals,[2],AllRadians),
	BestCorrCount = 1*length(AllRadians)*4/4,
	
	%ESetAllK4 = rotated_line_edges_coalesce(GrayVals,12,5, 4/4),
	%ESetAllK2 = rotated_line_edges_coalesce(GrayVals,12,5, 2/4),
	
	write_img(ImgOut,W,H,  
			  fun(P) -> 
					  case dict:find(P,EDict) of 
						  error -> 255; 
						  {ok,Value} -> round(255 - 255*erlang:min(1, Value/BestCorrCount))  
					  end 
			  end).
			  %fun(P) -> set_element_map(P,255,[{ESetAllK4,0},{ESetAllK2,128}]) end).

rotated_line_edges_coalesce_sum(GrayVals,AvgVals,AllRadians) ->
	lists:foldl(fun(V,EDict) ->
						EDictV = rotated_line_edges_coalesce(GrayVals,V,AllRadians),
						dict_sum(EDict,EDictV)
				end, dict:new(), AvgVals).

set_element_map(_P,DefVal,[]) -> DefVal;
set_element_map(P,DefVal,[{Set,Val}|Tail]) -> 
	case sets:is_element(P,Set) of true -> Val; false -> set_element_map(P,DefVal,Tail) end.

all_radians(DeltaLine) ->
    [ G / 180 * math:pi() || G <- lists:seq(0,179,DeltaLine) ].

rotated_line_edges_coalesce(GraysVals,AvgSize,AllRadians) ->
	%DirCount = length(AllRadians),
	DirProcess = fun(LinesRadians,Edges) ->
		io:format("Computing line radians ~p~n",[LinesRadians]),
		EOrig = all_rotated_line_edges(GraysVals,AvgSize,LinesRadians),
		RoundList = [ P || {X,Y} <- EOrig, P <- int_points_around({X,Y}) ],
		dict_recount(Edges,RoundList)								 
	end,
	_Edges = lists:foldl(DirProcess, dict:new(), AllRadians).
%	dict:fold(fun(Point,Count,ESet) ->
%					  case Count > DirCount*MinKCorrelate of
%						  true -> sets:add_element(Point, ESet);
%						  false -> ESet
%					  end
%	end,  sets:new(), Edges).

dict_recount(Dict,NewVals) ->
	lists:foldl(fun(V,D) -> dict:update_counter(V, 1, D) end, Dict, NewVals).

dict_sum(Dict1,Dict2) ->
	dict:fold(fun(V,Count,D) -> dict:update_counter(V, Count, D) end, Dict1, Dict2).

int_points_around({X,Y}) ->
	X0=trunc(X),
	Y0=trunc(Y),
	[ {X0,Y0},{X0+1,Y0},{X0,Y0+1},{X0+1,Y0+1} ].

line_edges(GrayVals,Y0,AvgSize) ->
	Xs = lists:seq(AvgSize, GrayVals#gray_vals.w - 1 - AvgSize),
	Vals = [ {X,avg_horiz_diff(GrayVals,X,Y0,AvgSize)} || X <- Xs ],
	%lists:foreach(fun({X,A}) ->
	%    io:format("~p  ~p~n",[X,trunc(A)])
	%end, Vals),
	dominant_vals(Vals,AvgSize).

lists_seq_safe(First,Last) ->
	case Last >= First of
		true -> lists:seq(First,Last);
		false -> []
	end.

dominant_vals(Vals,DominateSize) when is_list(Vals) ->
	Tup = list_to_tuple(Vals),
	IDoms = lists:filter( fun(I) ->
								  {_Locat1,AI} = element(I,Tup),
								  NeighTup = [ element(J,Tup) || J <- lists:seq(I-DominateSize, I+DominateSize) -- [I] ],
								  Neigh = [ A || {_Locat,A} <- NeighTup ],
								  AI > lists:max(Neigh)
						  end, lists_seq_safe(DominateSize+1, size(Tup)-DominateSize) ),
	[ element(I,Tup)  || I <- IDoms ].


avg_horiz_diff(GrayVals,X0,Y0,AvgSize) ->
	Left = lists:sum([ pixel(GrayVals,X,Y0) || X <- lists:seq(X0-AvgSize,X0-1) ])/AvgSize,
	Right = lists:sum([ pixel(GrayVals,X,Y0) || X <- lists:seq(X0+1,X0+AvgSize) ])/AvgSize,
	abs(Left-Right).



all_rotated_line_edges(GrayVals,AvgSize,LinesRadians) ->
	Center = {GrayVals#gray_vals.w/2, GrayVals#gray_vals.h/2},
	AcrossRadians = LinesRadians + math:pi()/2,
	DirAcross = {math:cos(AcrossRadians), -math:sin(AcrossRadians)},
	DirAlong = {math:cos(LinesRadians), -math:sin(LinesRadians)},
	Max = GrayVals#gray_vals.w + GrayVals#gray_vals.h,
	LineStarts = [ pshift(Center,DirAcross,I) || I <- lists:seq(-Max,Max) ],
	[ XY || P0 <- LineStarts, {XY,_Val} <- rotated_line_edges(GrayVals,AvgSize,P0, DirAlong) ].


rotated_line_edges(GrayVals,AvgSize,{X0,Y0}=_XY0, {Dx,Dy}=_DirXY) ->
	Max = GrayVals#gray_vals.w + GrayVals#gray_vals.h,
	%io:format("mpoints~w~n",[{X0 + Max*Dx, Y0 + I*Dy}]),
	XYAll = [ {X0 + I*Dx, Y0 + I*Dy} || I <- lists:seq(-Max, Max), in_bounds(GrayVals,X0 + I*Dx, Y0 + I*Dy) ],
	%    io:format("xy points~w~n",[XYAll]),
	%    [ {io:format("~w~n",[{X,Y}]), io:format("~w~n",[pixel(GrayVals,X,Y)])} || {X,Y} <- XYAll ],
	Points = [ {{X,Y},pixel(GrayVals,X,Y)} || {X,Y} <- XYAll ],
	%    io:format("points~w~n",[Points]),
	PointsTup = list_to_tuple(Points),
	IFirst = AvgSize+1,
	ILast = size(PointsTup)-AvgSize,
	case ILast >= IFirst of
		true ->
			PointsVals = [ {XY,avg_diff_along_line(PointsTup,I,AvgSize)} ||
						   I <- lists:seq(IFirst, ILast), {XY,_Val} <- [element(I,PointsTup)] ],
			dominant_vals(PointsVals, AvgSize);
		false -> [] % nothing to analyze
	end.


avg_diff_along_line(PointsTup,Idx,AvgSize) ->
	Left  = lists:sum([ Val || I <- lists:seq(Idx-AvgSize,Idx-1), {_XY,Val} <- [element(I,PointsTup)] ])/AvgSize,
	Right = lists:sum([ Val || I <- lists:seq(Idx+1,Idx+AvgSize), {_XY,Val} <- [element(I,PointsTup)] ])/AvgSize,
	abs(Left-Right).

%%% =========================================  Pixel array  =========================================

pshift(_Point={X,Y},_Dir={Dx,Dy},Times) ->
	{X+Times*Dx, Y+Times*Dy}.

in_bounds(GrayVals,X,Y) ->
	X>=0 andalso Y>=0 andalso X =< GrayVals#gray_vals.w-1 andalso Y =< GrayVals#gray_vals.h-1.

pixel(GrayVals,X,Y) when is_integer(X), is_integer(Y) ->
	case in_bounds(GrayVals,X,Y) of
		true ->
			element(1 + X + Y * GrayVals#gray_vals.w, GrayVals#gray_vals.byte_tuple);
		false ->
			erlang:error({out_of_image_bounds,X,Y})
	end;

pixel(GrayVals,X,Y) ->
	{X1,Y1} = {trunc(X),trunc(Y)},
	{X2,Y2} = {X1+1, Y1+1},
	P11 = pixel(GrayVals,X1,Y1),
	P12 = case Y==Y1 of true -> 0; false -> pixel(GrayVals,X1,Y2) end, % get pixel only if neeeded
	P21 = case X==X1 of true -> 0; false -> pixel(GrayVals,X2,Y1) end, % get pixel only if neeeded
	P22 = case (X==X1 orelse Y==Y1) of true -> 0; false -> pixel(GrayVals,X2,Y2) end, % get pixel only if neeeded
	PY1 = P21 * (X-X1) + P11 * (X2 - X),
	PY2 = P22 * (X-X1) + P12 * (X2 - X),
	PY2 * (Y-Y1) + PY1 * (Y2-Y).

%%% =========================================  I/O  =========================================

log_block_subpixel(GrayVals,X0,Y0) ->
	lists:foreach(fun(Y) ->
						  lists:foreach(fun(X) ->
												G = pixel(GrayVals,X0 + X/10, Y0 + Y/10 ),
												io:format("~4w",[ round(G) ])
										end, lists:seq(0,10)),
						  io:format("~n")
				  end, lists:seq(0,10)).

log_block(GrayVals,X0,Y0,Size) ->
	lists:foreach(fun(Y) ->
						  lists:foreach(fun(X) ->
												G = pixel(GrayVals,X,Y),
												io:format("~4w",[G])
										end, lists:seq(X0-Size,X0+Size)),
						  io:format("~n")
				  end, lists:seq(Y0-Size,Y0+Size)).

all_rotated_line_edges_gray_fun(GraysVals,AvgSize,LinesRadians) ->
	RoundList = [ {round(X),round(Y)} || {X,Y} <- all_rotated_line_edges(GraysVals,AvgSize,LinesRadians) ],
	ESet = sets:from_list(RoundList),
	fun(P) -> case sets:is_element(P,ESet) of true -> 0; false -> 255 end end.

write_img(File,W,H,GrayFun) when is_list(File) ->
	Xs = lists:seq(0,W-1),
	Port = open_port({spawn,"convert -size " ++ integer_to_list(W) ++ "x" ++ integer_to_list(H)
						  ++ " -depth 8 gray:- " ++ File},[]),
	lists:foreach(fun(Y) ->
						  Xs = lists:seq(0,W-1),
						  Line = << << V >> || X <- Xs, V <- [GrayFun({X,Y})] >>,
						  Port ! {self(), {command, Line}}
				  end, lists:seq(0,H-1)),
	port_close(Port).

write_img_edges(GrayVals,AvgSize) ->
	W = GrayVals#gray_vals.w,
	H = GrayVals#gray_vals.h,
	Xs = lists:seq(0,W-1),
	Port = open_port({spawn,"convert -size " ++ integer_to_list(W) ++ "x" ++ integer_to_list(H)
						  ++ " -depth 8 gray:- medge.png"},[]),
	lists:foreach(fun(Y) ->
						  
						  ESet = sets:from_list([ X || {X,_Diff} <- line_edges(GrayVals,Y,AvgSize) ]),
						  Xs = lists:seq(0,W-1),
						  Line = << << case sets:is_element(X,ESet) of true -> 0; false -> 255 end >> || X <- Xs >>,
Port ! {self(), {command, Line}}

end, lists:seq(0,H-1)),

port_close(Port).

write_img_lines(GrayVals,Y) ->
W = GrayVals#gray_vals.w,
H = 100,
Xs = lists:seq(0,W-1),
Line = << <<(pixel(GrayVals,X,Y))>> || X <- Xs >>,
Port = open_port({spawn,"convert -size " ++ integer_to_list(W) ++ "x" ++ integer_to_list(H)
++ " -depth 8 gray:- mlines.png"},[]),
[ Port ! {self(), {command, Line}} || _I <- lists:seq(1,H) ],
port_close(Port).

write_img_lines_example() ->
Port = open_port({spawn,"convert -size 3x2 -depth 8 gray:- m1.png"},[]),
Port ! {self(), {command, <<200,255,255,100,0,0>>}},
port_close(Port).

port_recv_all(Cmd) ->
Port = open_port({spawn,Cmd},[exit_status]),
lists:flatten( port_recv_all(Port,[]) ).

port_recv_all(Port,Acc) ->
receive
{Port,{exit_status,Status}} ->
io:format("Exit status ~p~n",[Status]),
lists:reverse(Acc);
{Port,{data,X}} ->
port_recv_all(Port,[X|Acc])
end.
